##########################################################################################

setwd("~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime")

##########################################################################################

library(ggplot2)
library(tidyverse)
library(readxl)
library(magrittr)

##########################################################################################

input_file <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime/ipmn-timing-master/data/Timing_Metrics.xlsx"
out_path <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime"
code_path <- "~/20220915_gastric_multiple/dna_combinePublic/scripts/evolutionTime/ipmn-timing-master"

##########################################################################################

source(paste0(code_path, "/code/functions.R"))

##########################################################################################

unixFriendly <- function(x){
  x <- colnames(x)
  x <- tolower(x)
  x <- gsub(" ", "_", x)
  x
}

##########################################################################################

dat <- read_excel(input_file,
                  sheet=2) %>%
  "["(-nrow(.), ) %>%
  set_colnames(unixFriendly(.)) %>%
  mutate(additional=mutations_ca-mutations_hg_ipmn)
avgN <- dat %>%
  summarize(mean=mean(additional)) %>%
  "[["(1)


##########################################################################################

ipmn.age <- read_excel(input_file,
                   sheet=1) %>%
"["(1:2, ) %>%
set_colnames(c("age", "n", "url")) %>%
mutate(age=as.numeric(age),
     n=as.integer(n))
pdac.age <- read_excel(input_file,
                   sheet=1) %>%
"["(5:6, ) %>%
set_colnames(c("age", "n", "url")) %>%
mutate(age=as.numeric(age),
     n=as.integer(n))
Tavg <- sum((pdac.age$age * pdac.age$n))/sum(pdac.age$n) -
sum(ipmn.age$age * ipmn.age$n)/sum(ipmn.age$n)
round(Tavg, 2)
mu <- avgN/Tavg ## mutations per year

##########################################################################################
## Prior model 1 (middle)
## mu*Tavg ## 15.5
## posterior is also gamma.  center mean at a/b=15.5
## The shape parameter $a$ can be interpreted as the sum of counts from b prior observations.  
## Below, we assume `r a=7.75; a` mutations in `r b=0.5; b` observations. 
## This prior is centered at `r round(a/b, 1)` mutations.
b <- 0.5
a <- avgN*b

aprime <- a + dat$additional
bprime <- b + 1  ## one sample
posterior <- matrix(NA, 1000, nrow(dat))
for(j in seq_along(aprime)){
  posterior[, j] <- rgamma(1000, aprime[j], rate=bprime)
}

avg.years <- colMeans(posterior) / mu
lower <- apply(posterior, 2, quantile, 0.025) / mu
upper <- apply(posterior, 2, quantile, 0.975) / mu

##########################################################################################
## Prior model 2
## mu*Tavg ## 15.5
## posterior is also gamma.  center mean at a/b=15.5
b <- 1
a <- avgN*b

aprime <- a + dat$additional
bprime <- b + 1  ## one sample
posterior <- matrix(NA, 1000, nrow(dat))
for(j in seq_along(aprime)){
  posterior[, j] <- rgamma(1000, aprime[j], rate=bprime)
}
avg.years2 <- colMeans(posterior) / mu
lower2 <- apply(posterior, 2, quantile, 0.025) / mu
upper2 <- apply(posterior, 2, quantile, 0.975) / mu

##########################################################################################
## Prior model 3 (longest IPMN->PDAC)
b <- 1/4
a <- avgN*b

aprime <- a + dat$additional
bprime <- b + 1  ## one sample
posterior <- matrix(NA, 1000, nrow(dat))
for(j in seq_along(aprime)){
  posterior[, j] <- rgamma(1000, aprime[j], rate=bprime)
}
avg.years3 <- colMeans(posterior) / mu
lower3 <- apply(posterior, 2, quantile, 0.025) / mu
upper3 <- apply(posterior, 2, quantile, 0.975) / mu

##########################################################################################
## Summary

priormodels <- rep(paste0("prior_model_", 1:3), each=nrow(dat)) %>%
  factor(., levels=unique(.))
facetlabels <- gsub("p", "P", levels(priormodels)) %>%
  gsub("_", " ", .)
facetlabels <- setNames(facetlabels, levels(priormodels))
dat$case <- factor(dat$case, levels=rev(dat$case))
tab <- tibble(id=rep(dat$case, 3),
              years=c(avg.years, avg.years2, avg.years3),
              lower=c(lower, lower2, lower3),
              upper=c(upper, upper2, upper3)) %>%
  mutate(id=factor(id, levels=rev(dat$case))) %>%
  mutate(a=rep(7.75, nrow(.)),
         b=rep(c(0.5, 1, 0.025), each=nrow(dat)),
         prior_model=priormodels) %>%
  left_join(dat, by=c("id"="case"))
p <- ggplot(tab, aes(ymin=lower, ymax=upper, x=id, xend=id)) +
  geom_point(aes(y=years), shape=21, color="gray") +
  geom_errorbar(width=0.2) +
  coord_flip() +
  facet_wrap(~prior_model, labeller=labeller(prior_model=facetlabels)) +
  theme(panel.background=element_rect(fill="white", color="black")) +
  xlab("") + ylab("Years")

out_name <- paste0(out_path , "/timing.pdf")
ggsave( out_name , p)


## Averages
targeted <- tibble(rev(levels(tab$id)),
                   captured=c(1,1,1,0,1,
                              0,0,0,1,0,
                              0,0,0,1,0,
                              1,0,0)) %>%
  set_colnames(c("id", "captured")) %>%
  mutate(id=factor(id, levels=levels(tab$id)))
tab2 <- left_join(tab, targeted, by="id")
tab2 %>% group_by(prior_model) %>%
  summarize(average=mean(years),
            lower=mean(lower),
            upper=mean(upper))

tab2 %>% group_by(prior_model, captured) %>%
  summarize(average=mean(years),
            lower=mean(lower),
            upper=mean(upper)) %>%
  arrange(captured)